Microbiome assessment for A. cervicornis microplastics project (16S and ITS2)
Last run on 02 August 2022
16S Assessment
I am largely following the DADA2 tutorial by Benjamin Callahan that can be found here.
Inspecting the data
Check samples to use
## [1] "A10" "A11" "A12" "A2" "A3" "A6" "A8" "A9" "B10" "B11" "B3" "B4"
## [13] "B5" "B7" "C1" "C10" "C11" "C2" "C4" "G5" "G6"
View example quality plots
Viewing the quality of 2 of the forward samples 
Viewing the quality of the same 2 of the reverse samples 
Filter and trim samples
Perform the filtering
You can see in above figures that the quality of reads drop off towards the end, so we need to filter out these low quality reads
## reads.in reads.out
## A10_16S_R1.fastq 25343 23991
## A11_16S_R1.fastq 27398 26006
## A12_16S_R1.fastq 45822 44152
## A2_16S_R1.fastq 99390 94836
## A3_16S_R1.fastq 77568 73174
## A6_16S_R1.fastq 44719 41897After filtering out the low quality reads, we have maintained about 95.7 of the original reads.
Check the filtering
Viewing the quality of the same 2 of the forward samples post filtering and trimming 
Viewing the quality of the same 2 of the reverse samples post filtering and trimming 
DADA2 sample processing
Error Rates
The error rates for each possible transition (A→C, A→G, …) are shown. Points are the observed error rates for each consensus quality score. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. We want the estimated error rates (black line) to be a good fit to the observed rates (points), and the error rates to drop with increased quality.
Sample Inference
We are now ready to apply the core sample inference algorithm to the filtered and trimmed sequence data.
Forward sample inference
## Sample 1 - 23991 reads in 6094 unique sequences.
## Sample 2 - 26006 reads in 6312 unique sequences.
## Sample 3 - 44152 reads in 7218 unique sequences.
## Sample 4 - 94836 reads in 13509 unique sequences.
## Sample 5 - 73174 reads in 9644 unique sequences.
## Sample 6 - 41897 reads in 7509 unique sequences.
## Sample 7 - 48504 reads in 7342 unique sequences.
## Sample 8 - 32048 reads in 7630 unique sequences.
## Sample 9 - 24910 reads in 5725 unique sequences.
## Sample 10 - 34782 reads in 5737 unique sequences.
## Sample 11 - 34374 reads in 5211 unique sequences.
## Sample 12 - 90460 reads in 8319 unique sequences.
## Sample 13 - 23566 reads in 5815 unique sequences.
## Sample 14 - 110660 reads in 10438 unique sequences.
## Sample 15 - 27927 reads in 5388 unique sequences.
## Sample 16 - 73356 reads in 7647 unique sequences.
## Sample 17 - 23594 reads in 6498 unique sequences.
## Sample 18 - 25590 reads in 6768 unique sequences.
## Sample 19 - 42054 reads in 7901 unique sequences.
## Sample 20 - 1795 reads in 421 unique sequences.
## Sample 21 - 162 reads in 63 unique sequences.Inspecting the returned dada-class object for the first forward sample:
## dada-class: object describing DADA2 denoising results
## 113 sequence variants were inferred from 6094 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16The DADA2 algorithm inferred 144 true sequence variants from the 9515 unique sequences in the first sample. There is much more to the dada-class return object than this (see help("dada-class") for some info), including multiple diagnostics about the quality of each denoised sequence variant, but that is beyond the scope of an introductory tutorial.
Reverse sample inference
## Sample 1 - 23991 reads in 5924 unique sequences.
## Sample 2 - 26006 reads in 6092 unique sequences.
## Sample 3 - 44152 reads in 7108 unique sequences.
## Sample 4 - 94836 reads in 12281 unique sequences.
## Sample 5 - 73174 reads in 8686 unique sequences.
## Sample 6 - 41897 reads in 7149 unique sequences.
## Sample 7 - 48504 reads in 7069 unique sequences.
## Sample 8 - 32048 reads in 7283 unique sequences.
## Sample 9 - 24910 reads in 5526 unique sequences.
## Sample 10 - 34782 reads in 5706 unique sequences.
## Sample 11 - 34374 reads in 4744 unique sequences.
## Sample 12 - 90460 reads in 8211 unique sequences.
## Sample 13 - 23566 reads in 5706 unique sequences.
## Sample 14 - 110660 reads in 10194 unique sequences.
## Sample 15 - 27927 reads in 5233 unique sequences.
## Sample 16 - 73356 reads in 7528 unique sequences.
## Sample 17 - 23594 reads in 6469 unique sequences.
## Sample 18 - 25590 reads in 6315 unique sequences.
## Sample 19 - 42054 reads in 7772 unique sequences.
## Sample 20 - 1795 reads in 465 unique sequences.
## Sample 21 - 162 reads in 69 unique sequences.
Merge paired reads
We now merge the forward and reverse reads together to obtain the full denoised sequences. Merging is performed by aligning the denoised forward reads with the reverse-complement of the corresponding denoised reverse reads, and then constructing the merged “contig” sequences. By default, merged sequences are only output if the forward and reverse reads overlap by at least 12 bases, and are identical to each other in the overlap region (but these conditions can be changed via function arguments).
The mergers object is a list of data.frames from each sample. Each data.frame contains the merged sequence, its abundance, and the indices of the forward and reverse sequence variants that were merged. Paired reads that did not exactly overlap were removed by mergePairs, further reducing spurious output.
Construct sequence table
We can now construct an amplicon sequence variant table (ASV) table, a higher-resolution version of the OTU table produced by traditional methods.
##
## 223 227 239 246 251 252 253 254 255 256 257
## 2 1 1 1 5 46 1601 86 4 4 2
The sequence table is a matrix with rows corresponding to (and named by) the samples, and columns corresponding to (and named by) the sequence variants. This table contains 1753 ASVs.
After viewing the distribution of read lengths, it looks like we have some that fall outside the expected range (244 - 264) so we will go ahead and remove these non target length sequences.
##
## 246 251 252 253 254 255 256 257
## 1 5 46 1601 86 4 4 2This updated table now contains 1749 ASVs.
Remove chimeras
The core dada method corrects substitution and indel errors, but chimeras remain. Fortunately, the accuracy of sequence variants after denoising makes identifying chimeric ASVs simpler than when dealing with fuzzy OTUs. Chimeric sequences are identified if they can be exactly reconstructed by combining a left-segment and a right-segment from two more abundant “parent” sequences.
A total of 9 bimeras were identified from the 1740 input sequences, thus retaining 100% of sequences.
Track reads through the pipeline
As a final check of our progress, we can look at the number of reads that made it through each step in the pipeline:
| Genotype | Treatment | Input | Filtered | Denoised Forward | Denoised Reverse | Merged | Nonchimera |
|---|---|---|---|---|---|---|---|
| 6 | OAW MP+ | 25343 | 23991 | 23882 | 23969 | 23811 | 23811 |
| 5 | OAW Control | 27398 | 26006 | 25982 | 25971 | 25961 | 25961 |
| 5a | AMB MP+ | 45822 | 44152 | 44079 | 44127 | 43982 | 43982 |
| 10a | OAW MP+ | 99390 | 94836 | 94654 | 94676 | 94353 | 94205 |
| 5a | AMB Control | 77568 | 73174 | 73129 | 73130 | 73018 | 73018 |
| 10a | AMB MP+ | 44719 | 41897 | 41808 | 41848 | 41496 | 41496 |
| 6 | AMB Control | 50868 | 48504 | 48363 | 48414 | 48114 | 48073 |
| 5a | OAW Control | 33629 | 32048 | 31998 | 32015 | 31874 | 31842 |
| 5a | OAW MP+ | 26111 | 24910 | 24868 | 24838 | 24666 | 24641 |
| 5 | AMB Control | 36048 | 34782 | 34765 | 34754 | 34710 | 34710 |
| 2 | AMB MP+ | 36469 | 34374 | 34349 | 34353 | 34289 | 34289 |
| 2 | AMB Control | 93851 | 90460 | 90421 | 90323 | 90201 | 90201 |
| 6 | AMB MP+ | 24767 | 23566 | 23506 | 23540 | 23364 | 23328 |
| 10a | OAW Control | 114866 | 110660 | 110570 | 110599 | 110465 | 110452 |
| 5 | AMB MP+ | 29051 | 27927 | 27900 | 27896 | 27851 | 27851 |
| 2 | OAW MP+ | 75311 | 73356 | 73246 | 73255 | 72852 | 72852 |
| 6 | OAW Control | 24571 | 23594 | 23583 | 23578 | 23571 | 23571 |
| 10a | AMB Control | 26776 | 25590 | 25551 | 25561 | 25491 | 25491 |
| 5 | OAW MP+ | 43598 | 42054 | 42021 | 42021 | 41807 | 41723 |
| NA | blank_control | 1918 | 1795 | 1792 | 1792 | 1792 | 1792 |
| NA | blank_control | 169 | 162 | 153 | 158 | 153 | 153 |
Assign taxonomy
It is common at this point, especially in 16S/18S/ITS amplicon sequencing, to assign taxonomy to the sequence variants. The DADA2 package provides a native implementation of the naive Bayesian classifier method for this purpose. The assignTaxonomy function takes as input a set of sequences to be classified and a training set of reference sequences with known taxonomy, and outputs taxonomic assignments with at least minBoot bootstrap confidence.
The dada2 package GitHub maintains the most updated versions of the Silva databases., but I downloaded the databases from the associated Zenodo. The versions in this GitHub repository, used here, were last updated on 26 July 2022.
Let’s inspect the taxonomic assignments:
## Kingdom Phylum Class Order
## [1,] "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rickettsiales"
## [2,] "Bacteria" "Cyanobacteria" "Cyanobacteriia" "Chloroplast"
## [3,] "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Pseudomonadales"
## [4,] "Bacteria" "Myxococcota" "Myxococcia" "Myxococcales"
## [5,] "Bacteria" "Campylobacterota" "Campylobacteria" "Campylobacterales"
## [6,] "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Burkholderiales"
## Family Genus Species
## [1,] "Fokiniaceae" "MD3-55" NA
## [2,] NA NA NA
## [3,] "Pseudomonadaceae" "Pseudomonas" NA
## [4,] "Myxococcaceae" "P3OB-42" NA
## [5,] NA NA NA
## [6,] "Alcaligenaceae" "Alcaligenes" NAGreat! We can now save this and hand it off to Phyloseq for further analyses.
Remove contaminations
Remove mitochondria, chloroplasts, and non-bacteria
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1610 taxa and 21 samples ]
## sample_data() Sample Data: [ 21 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 1610 taxa by 7 taxonomic ranks ]Removal of mitochondira, chloroplasts, and non-bacteria taxa reduced the total number of taxa from 1740 to 1610 (130 total).
Remove neg control contamination
Just 3 out of the 1610 ASVs were classified as contaminants.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1607 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 1607 taxa by 7 taxonomic ranks ]
Blast NCBI
I have not rerun this portion since only running the pipeline with the 19 GE samples!
Looks good! We have now cleaned up the sample data to remove contamination from non target organisms and those from the negative controls.
Rarefy
Rarefy cleaned data
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1538 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 1538 taxa by 8 taxonomic ranks ]## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1509 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 1509 taxa by 8 taxonomic ranks ]
Trim underrepresented OUTs
## [1] "samples with counts below z-score -2.5 :"
## character(0)
## [1] "zscores:"
## named numeric(0)
## [1] "OTUs passing frequency cutoff 1e-04 : 416"
## [1] "OTUs with counts in 0.02 of samples:"
##
## TRUE
## 416## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 416 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 416 taxa by 8 taxonomic ranks ]
Rarefy trimmed data
## named numeric(0)
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 416 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 416 taxa by 8 taxonomic ranks ]
Microbiome diversity analyses
Notes from phyloseq author Visualize alpha-diversity - Should be done on raw, untrimmed dataset
Simpson diversity
## OK: residuals appear as normally distributed (p = 0.496).## OK: Error variance appears to be homoscedastic (p = 0.428).## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: InvSimpson
## Chisq Df Pr(>Chisq)
## Treatment 5.7787 3 0.1229## Linear mixed model fit by REML ['lmerMod']
## Formula: InvSimpson ~ Treatment + (1 | Genotype)
## Data: alpha_df_div
##
## REML criterion at convergence: 66.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5662 -0.4630 -0.1005 0.4098 1.5485
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 1.404 1.185
## Residual 2.453 1.566
## Number of obs: 19, groups: Genotype, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.7669 0.8783 2.012
## TreatmentAMB MP+ 0.9121 0.9905 0.921
## TreatmentOAW Control 2.5412 1.0628 2.391
## TreatmentOAW MP+ 0.9365 0.9905 0.945
##
## Correlation of Fixed Effects:
## (Intr) TAMBMP TrOAWC
## TrtmnAMBMP+ -0.564
## TrtmntOAWCn -0.526 0.466
## TrtmnOAWMP+ -0.564 0.500 0.466Shannon diversity
## OK: residuals appear as normally distributed (p = 0.883).## OK: Error variance appears to be homoscedastic (p = 0.917).## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Shannon
## Chisq Df Pr(>Chisq)
## Treatment 4.2029 3 0.2404## Linear mixed model fit by REML ['lmerMod']
## Formula: Shannon ~ Treatment + (1 | Genotype)
## Data: alpha_df_div
##
## REML criterion at convergence: 47.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.63835 -0.64576 -0.09138 0.56195 1.58016
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 0.2279 0.4774
## Residual 0.7302 0.8545
## Number of obs: 19, groups: Genotype, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.1350 0.4377 2.593
## TreatmentAMB MP+ 0.7171 0.5404 1.327
## TreatmentOAW Control 1.1529 0.5783 1.994
## TreatmentOAW MP+ 0.5053 0.5404 0.935
##
## Correlation of Fixed Effects:
## (Intr) TAMBMP TrOAWC
## TrtmnAMBMP+ -0.617
## TrtmntOAWCn -0.577 0.467
## TrtmnOAWMP+ -0.617 0.500 0.467Species richness
## OK: residuals appear as normally distributed (p = 0.060).## OK: Error variance appears to be homoscedastic (p = 0.429).## Anova Table (Type II tests)
##
## Response: log(Observed)
## Sum Sq Df F value Pr(>F)
## Treatment 0.05346 3 0.2717 0.8449
## Residuals 0.98389 15##
## Call:
## lm(formula = log(Observed) ~ Treatment, data = alpha_df_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32518 -0.15924 -0.06924 0.14390 0.59111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.60184 0.11454 40.178 <2e-16 ***
## TreatmentAMB MP+ -0.02151 0.16198 -0.133 0.896
## TreatmentOAW Control 0.12179 0.17180 0.709 0.489
## TreatmentOAW MP+ 0.05102 0.16198 0.315 0.757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2561 on 15 degrees of freedom
## Multiple R-squared: 0.05153, Adjusted R-squared: -0.1382
## F-statistic: 0.2717 on 3 and 15 DF, p-value: 0.8449Evenness
## OK: residuals appear as normally distributed (p = 0.946).## OK: Error variance appears to be homoscedastic (p = 0.789).## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(Evenness)
## Chisq Df Pr(>Chisq)
## Treatment 3.759 3 0.2887## Linear mixed model fit by REML ['lmerMod']
## Formula: log(Evenness) ~ Treatment + (1 | Genotype)
## Data: alpha_df_div
##
## REML criterion at convergence: 38.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.40328 -0.65373 -0.00822 0.50024 1.79776
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 0.1702 0.4126
## Residual 0.4009 0.6332
## Number of obs: 19, groups: Genotype, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.6471 0.3380 -4.874
## TreatmentAMB MP+ 0.6716 0.4005 1.677
## TreatmentOAW Control 0.6933 0.4291 1.615
## TreatmentOAW MP+ 0.3548 0.4005 0.886
##
## Correlation of Fixed Effects:
## (Intr) TAMBMP TrOAWC
## TrtmnAMBMP+ -0.592
## TrtmntOAWCn -0.553 0.467
## TrtmnOAWMP+ -0.592 0.500 0.467All alpha diversity metrics
Core vs. accessory microbiome
Core
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 11 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 11 taxa by 8 taxonomic ranks ]
Core stats
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.02979 0.0099296 0.4555 0.7173
## Residuals 15 0.32697 0.0217981
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
##
## adonis2(formula = seq_core ~ Treatment + Genotype, data = samdf_core, permutations = 1500)
## Df SumOfSqs R2 F Pr(>F)
## Treatment 3 0.14373 0.17424 1.2359 0.3371
## Genotype 4 0.25476 0.30885 1.6431 0.2085
## Residual 11 0.42640 0.51691
## Total 18 0.82489 1.00000## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.010069 0.0033562 0.2065 0.8903
## Residuals 15 0.243831 0.0162554
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
##
## adonis2(formula = seq_core_rel ~ Treatment + Genotype, data = samdf_core, permutations = 1500)
## Df SumOfSqs R2 F Pr(>F)
## Treatment 3 0.07444 0.14086 0.9785 0.4177
## Genotype 4 0.17508 0.33129 1.7260 0.1972
## Residual 11 0.27896 0.52785
## Total 18 0.52848 1.00000
Accessory
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 405 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 6 sample variables ]
## tax_table() Taxonomy Table: [ 405 taxa by 8 taxonomic ranks ]
Accessory stats
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.0011202 0.00037341 0.3049 0.8214
## Residuals 15 0.0183684 0.00122456
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
##
## adonis2(formula = seq_acc ~ Treatment + Genotype, data = samdf_acc, permutations = 1500)
## Df SumOfSqs R2 F Pr(>F)
## Treatment 3 1.1413 0.14099 0.8319 0.9660
## Genotype 4 1.9230 0.23756 1.0512 0.2685
## Residual 11 5.0305 0.62145
## Total 18 8.0947 1.00000## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.0005439 0.00018130 0.1934 0.8992
## Residuals 15 0.0140613 0.00093742
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
##
## adonis2(formula = seq_acc_rel ~ Treatment + Genotype, data = samdf_acc, permutations = 1500)
## Df SumOfSqs R2 F Pr(>F)
## Treatment 3 1.0908 0.13698 0.8047 0.9820
## Genotype 4 1.9024 0.23890 1.0526 0.2898
## Residual 11 4.9701 0.62412
## Total 18 7.9633 1.00000
ITS2 Assessment
Abundance by treatment and genotype
Figure SXX. Relative abundance of major ITS2 types by coral fragment (x axis) and treatment (facets). Light green represent Symbiodinium spp. (A3) and dark green represents Breviolum spp. (B2).
Abundance by treatment only
Manuscript Materials
Methods
ITS2 and 16S metabarcoding
DNA was extracted using an RNAqueous kit (ThermoFisher Scientific) as described above, except samples were not subjected to the DNA removal step. Internal transcribed spacer region 2 (ITS2) PCR amplification was performed using SYM_VAR_5.8S2 and SYM_VAR_REV primers (as described in Hume et al. 2013, 2015) using the following PCR profile: 26 cycles of 95℃ for 40 s, 59℃ for 2 min, 72℃ for 1 min and a final extension of 72℃ for 7 min. A negative control was included and it failed to amplify so it was not sequenced. PCR products were cleaned using the GeneJET PCR Purification kit (ThermoFisher Scientific) according to the manufacturer’s instructions. A second PCR was performed to dual-barcode samples before pooling, which was done based on the visualization of band intensity on a 1% agarose gel. 20 μl of the pool was run on a 2% SYBR Green gel for 100 min at 70 V. The target band was excised and placed in 30 μl of Milli-Q water overnight at 4℃. The supernatant concentration was measured using a DeNovix DS-11+ Spectrophotometer.
The V4 region of the 16S rRNA gene was amplified via PCR using Hyb515f and Hyb806R primers (Parada et al. 2016, Apprill et al. 2015) and the following PCR profile: 30 cycles of 95℃ for 40 s, 63℃ for 2 min, 72℃ for 1 min and a final extension of 7 min. Two negative controls using Milli-Q water were also prepared during the 16S library preparations and these samples were also sequenced. The same procedure as described above for the ITS2 samples was then followed. Concentrations of the ITS2 and 16S pools were used to combine the two pools in a 1:3 ratio respectively. Libraries were sequenced on Illumina MiSeq (paired-end 250 bp) at Tufts Genomics Core Facility.
ITS2 and 16S analyses
Demultiplexed reads were pre-processed using bbmap (CITE) to remove adapters and to split ITS2 and 16S reads based on primers. Resulting ITS2 reads were then analyzed using the SymPortal framework (Hume et al., 2019) by submitting paired fastq.gz files directly to SymPortal. Samples were subjected to standard sequence quality control protocols using MOTHUR 1.39.5 (Schloss et al., 2009), the BLAST+ suite of executables (Camacho et al., 2009), and minimum entropy decomposition (Eren et al., 2015). SymPortal identifies specific sets of defining intragenomic ITS2 sequence variants (DIVs) to define ITS2 type profiles that are indicative of genetically differentiated Symbiodiniaceae taxa (Hume et al., 2019).
16S primers were removed using cutadapt (CITE) and DADA2 (CITE) was used to conduct quality filtering and inference of 1740 amplicon sequence variants (ASVs) (see Table SXX to track reads lost through filtering). 16S taxonomy was assigned using the Silva v. 138.1 database (CITE) and National Center for Biotechnology Information’s nucleotide database using blast+ (CITE). ASVs matching mitochondria, chloroplasts, or non-bacterial kingdoms were removed (130 total) and 3 ASVs were removed based on negative controls as contaminants (Decontam; (CITE)). Cleaned counts were then rarefied to 1.766810^{4} per sample using vegan (CITE), trimmed using MCM.OUT (CITE) to remove ASVs present in less than 0.01% of counts, and then rarefied again to 1.5610^{4} resulting in 416 ASVs. DIVERSITY AND COMMUNITY ANALYSES
Results
Results text goes here
Figures
Figure 2
Figure 2. Relative abundance of major ITS2 types grouped by treatment (n = 4-5 corals per treatment). Light purple represents Symbiodinium fitti (A3) and dark purple represents Breviolum minutum (B2).
Figure 5
Figure 5: Acropora cervicornis bacterial (16S) communities across experimental treatments including, (A) the core microbial taxa in each fragment, (B) Shannon diversity across treatments, and multivariate ordination plots (PCoA) of between-sample Bray–Curtis dissimilarity of rarefied ASVs of both (C) core and (D) accessory taxa. PCoA ellipses depict 95% confidence intervals and p-values indicate significance of dispersion (Pdisp) and multivariate location (Ptreat and Pgenotype).
Figure 6
Supplemental figures
Supplemental tables
Session Information
All code was written by Colleen B. Bove, feel free to contact with questions.
Session information from the last run date on 02 August 2022:
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] nlme_3.1-155 compositions_2.0-4
## [3] patchwork_1.1.1 lme4_1.1-28
## [5] car_3.1-0 carData_3.0-5
## [7] performance_0.8.0 awtools_0.2.1
## [9] microbiome_1.8.0 ggpubr_0.4.0
## [11] MCMC.OTU_1.0.10 MCMCglmm_2.33
## [13] ape_5.6-1 coda_0.19-4
## [15] Matrix_1.3-4 janitor_2.1.0
## [17] vegan_2.5-7 lattice_0.20-45
## [19] permute_0.9-7 plotly_4.10.0
## [21] RColorBrewer_1.1-3 decontam_1.6.0
## [23] phyloseq_1.30.0 kableExtra_1.3.4
## [25] ShortRead_1.44.3 GenomicAlignments_1.22.1
## [27] SummarizedExperiment_1.16.1 DelayedArray_0.12.3
## [29] matrixStats_0.62.0 Biobase_2.46.0
## [31] Rsamtools_2.2.3 GenomicRanges_1.38.0
## [33] GenomeInfoDb_1.22.1 Biostrings_2.54.0
## [35] XVector_0.26.0 IRanges_2.20.2
## [37] S4Vectors_0.24.4 BiocParallel_1.20.1
## [39] BiocGenerics_0.32.0 forcats_0.5.1
## [41] stringr_1.4.0 dplyr_1.0.8
## [43] purrr_0.3.4 readr_2.1.2
## [45] tidyr_1.2.0 tibble_3.1.7
## [47] ggplot2_3.3.6 tidyverse_1.3.1
## [49] dada2_1.20.0 Rcpp_1.0.9
## [51] knitr_1.33
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.4.1 systemfonts_1.0.2
## [4] plyr_1.8.7 igraph_1.2.11 lazyeval_0.2.2
## [7] splines_3.6.3 digest_0.6.29 foreach_1.5.2
## [10] htmltools_0.5.2 fansi_1.0.3 magrittr_2.0.3
## [13] cluster_2.1.2 tzdb_0.2.0 modelr_0.1.8
## [16] bayesm_3.1-4 RcppParallel_5.1.5 svglite_2.1.0
## [19] jpeg_0.1-9 colorspace_2.0-3 rvest_1.0.1
## [22] textshaping_0.3.6 haven_2.4.3 xfun_0.29
## [25] crayon_1.5.1 RCurl_1.98-1.7 jsonlite_1.7.2
## [28] survival_3.2-13 iterators_1.0.14 glue_1.6.2
## [31] gtable_0.3.0 zlibbioc_1.32.0 webshot_0.5.2
## [34] Rhdf5lib_1.8.0 DEoptimR_1.0-11 abind_1.4-5
## [37] scales_1.2.0 DBI_1.1.3 rstatix_0.7.0
## [40] viridisLite_0.4.0 htmlwidgets_1.5.4 httr_1.4.3
## [43] ellipsis_0.3.2 farver_2.1.1 pkgconfig_2.0.3
## [46] sass_0.4.0 dbplyr_2.1.1 deldir_1.0-6
## [49] utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1
## [52] rlang_1.0.4 reshape2_1.4.4 munsell_0.5.0
## [55] cellranger_1.1.0 tools_3.6.3 cachem_1.0.6
## [58] cli_3.3.0 generics_0.1.3 ade4_1.7-18
## [61] broom_1.0.0 evaluate_0.15 biomformat_1.14.0
## [64] fastmap_1.1.0 ragg_1.1.3 yaml_2.3.5
## [67] fs_1.5.2 robustbase_0.93-9 xml2_1.3.3
## [70] compiler_3.6.3 rstudioapi_0.13 png_0.1-7
## [73] ggsignif_0.6.3 reprex_2.0.1 bslib_0.4.0
## [76] stringi_1.7.8 highr_0.9 cubature_2.0.4.2
## [79] nloptr_1.2.2.3 tensorA_0.36.2 multtest_2.42.0
## [82] vctrs_0.4.1 pillar_1.8.0 lifecycle_1.0.1
## [85] jquerylib_0.1.4 cowplot_1.1.1 insight_0.18.0
## [88] data.table_1.14.2 bitops_1.0-7 corpcor_1.6.10
## [91] R6_2.5.1 latticeExtra_0.6-30 hwriter_1.3.2.1
## [94] gridExtra_2.3 codetools_0.2-18 boot_1.3-28
## [97] MASS_7.3-55 assertthat_0.2.1 rhdf5_2.30.1
## [100] withr_2.5.0 GenomeInfoDbData_1.2.2 mgcv_1.8-38
## [103] hms_1.1.1 grid_3.6.3 minqa_1.2.4
## [106] rmarkdown_2.13 snakecase_0.11.0 Rtsne_0.15
## [109] lubridate_1.8.0 interp_1.0-33